Disorder-Induced Stabilization of the Pseudogap in Strongly Correlated Systems 
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The interplay of strong interaction and strong disorder, as contained in the Anderson-Hubbard 
model, is addressed using two non-perturbative numerical methods: the Lanczos algorithm in the 
grand canonical ensemble at zero temperature and Quantum Monte Carlo. We find distinctive 
evidence for a zero-energy anomaly which is robust upon variation of doping, disorder and interaction 
strength. Its similarities to, and differences from, pseudogap formation in other contexts, including 
perturbative treatments of interactions and disorder, classical theories of localized charges, and in 
the clean Hubbard model, are discussed. 



"Pseudogap" anomalies in the single particle density 
of states are a central feature of seemingly disparate ma- 
terials and models. On the one hand, they have an early 
history in the metal-insulator transition and the study 
of the interplay between disorder and inter-particle in- 
teraction. In the metallic limit of weak coupling and 
weak disorder, Altshuler and Aronov (AA) showed[l(, by 
means of perturbation theory, that there is a depression 
of the spectral density at the chemical potential, the mag- 
nitude of the depression being dependent on the interac- 
tion strength. In the opposite limit of completely local- 
ized charges (where the model becomes classical), Efros 
and Shklovskii (ES)0 have shown that the combined ef- 
fect of the unscreened Coulomb potential and disorder 
also gives rise to an anomaly at the chemical potential - 
the Coulomb gap. 

On the other hand, in contrast to this situation in 
which randomness and electron-electron correlation both 
are crucial, pseudogap anomalies also arise in cupratc 
superconductors and the Hubbard model with no disor- 
der. While particle density does not play a central role 
in either the AA and ES pseudogaps, in the high T c ma- 
terials the pseudogap is confined to a low-doping region 
between the superconducting dome and the parent an- 
tiferromagnetic material. Likewise, in numerical studies 
of the Hubbard Hamiltonian 0, U, [B[ the pseudogap is 
absent for particle densities p < 0.80. 

That disorder is essential to the pseudogap in one sit- 
uation, yet present in the clean system in another, raises 
a fundamental question: What role does randomness 
play in low energy anomalies in the density of states of 
strongly correlated systems? In this manuscript, we sug- 
gest that randomness stabilizes the pseudogap. Indeed, 
we demonstrate two remarkable features of the pseudo- 
gap in the disordered Hubbard model. First, the density 
of states anomaly persists in the limit of an infinitely 
repulsive local potential U, even though the magnetic 
energy scale J oc t 2 /U is driven to zero. Second, it is 
independent of doping for a wide range of disorder and 
interaction strengths. The insensitivity of the pseudo- 
gap to doping is a novel effect related to the presence of 
disorder that has been observed experimentally^. 



We consider the Anderson-Hubbard Hamiltonian, 
H({ei}) = -tJ2'cl c i* + J2 6ini ° + \ U \J2 ni 1 ni l> W 



where the on-site energies e, are sampled uniformly from 
the interval [— A/2 : A/2], c\ a (ci a ) are fermion cre- 
ation(annihilation) operators for site i and spin a, and 
Uic = c\ a Ci a . The primed summation is on nearest neigh- 
bors only. We diagonalize this Hamiltonian on 10-site 
square clusters using the Lanczos algorithm at T = 
and on larger, 64-site clusters using finite temperature 
determinant Quantum Monte Carlo (DQMC). As usual, 
physical properties are computed as averages over many 
disorder realizations. 

Since both computational methods are well described 
in a number of previous publications 0, HI we focus 
our technical discussion only on the use of the grand- 
canonical ensemble to carry out the averaging process. 
In the DQMC case this is accomplished naturally as 
this method works, by construction, directly within a 
grand-canonical scheme. The Lanczos method, on the 
other hand, is a canonical technique that diagonalizcs the 
Hamiltonian in sectors of the Hilbcrt space with constant 
number of particles. Since we are interested in ground 
state properties, a given choice of chemical potential fi 
has the effect of singling out that particle number sector 
whose ground state minimizes (H — /j,N) . Which particu- 
lar sector is selected is, of course, dependent on the disor- 
der realization defining H ({ej}) and it is not known a pri- 
ori. Although this requires diagonalization of all particle 
sectors it also allows the treatment of non-commensurate 
fillings. Since the size of the system is fairly small we 
take advantage of the possibility of changing the bound- 
ary conditions of the electronic wave function to reduce 
finite size errors. The second quantized Hamiltonian be- 
comes therefore a function of k, a vector belonging to 
the first Brillouin zone of the simulation lattice. Within 
the grand-canonical ensemble k is treated as a quenched 
disorder variable and sampled uniformly. This is found 
to be important in the calculation of the spectral density 
for the model considered here. In the limit of vanish- 
ing disorder this scheme reduces to the integration over 
boundary condition tcchmque[9l. ficj| . With no disorder 
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FIG. 1: Evolution of the local spectral density, A(ui), at p — 1 
for U = 8 as a function of A. Left panel: A < U. The 
Mott gap evident for A = and 4 is partially filled in as the 
randomness increases, so that the value at the Fermi surface, 
A(lj = 0), becomes finite at A = U. Right panel: A > U. A 
pseudogap survives even for large randomness. The spectra 
are the average of 1000 independent disorder realizations. 

and no interaction one recovers the exact spectral density 
in the thermodynamic limit. 

We argue that the density of states obtained on the 
cluster sizes used here is relevant both because of the 
restoration of the thermodynamic limit by boundary con- 
dition averaging, and because finite size errors become in- 
creasingly small in the strong disorder case on which we 
focus. It is also worth noting that, unlike studies of other 
phenomena like long range magnetic or superconducting 
order in the Hubbard Hamiltonian, which critically rely 
on finite size scaling on large lattices, the density of states 
is less sensitive to system size. The AA calculation pro- 
vides a particularly clear example of this statement as 
the zero-energy anomaly was predicted using a class of 
diagrams that failed to capture localization. 

The disorder-averaged spectral density is given by 

i^i r 

A(Lo) = —Y,Ml^l dkd£MK{eihLo) (2) 

i— 1 

where N s and Qbz are, respectively, the number of sites 
and the volume of the Brillouin zone and the fermion 
addition part of Ai(k, {e^}, uj) is defined as 

MkAe i }M = - 1 -^(c i u; _ H l Eo+iri 4) (3) 

and computed using the continued fraction algorithm of 
Haydock et aZ.[H|]. The electron subtraction spectrum 
is obtained analogously by interchanging Cj and cj. In 
Eq. [3] Eq is the ground state energy of H = H(k, {£%}) 
and rj is a small real parameter giving the broadening 



FIG. 2: Evolution of the local spectral density away from 
half-filling (no Mott gap at A = 0). Left panel: A(ui) as a 
function of A. A ranges from well below U to well above. 
The dip in A{uj) at u) = becomes deeper as A increases. In 
the clean limit A — we observe no pseudogap this far from 
half-filling, in agreement with 

HHil]. Right panel: DQMC 
results on a 8x8=64 site lattice at T — t/6 averaging over 64 
disorder realizations. The scale (width) of the pseudogap is 
largely independent of doping. For p — 0.6 the 10-site Lanczos 
result is also reported. Despite the different cluster size the 
agreement between the two techniques is excellent. 

of the 8 functions constituting the spectrum. A(lo) is 
directly measured in tunneling spectroscopy, photoemis- 
sion and inverse photoemission experiments. Since disor- 
der averaging restores particle-hole symmetry, it suffices 
to consider the evolution of A(oj) at densities p < 1. In 
DQMC, A(w) is obtained by a maximum entropy anal ytic 
continuation of the imaginary time Green's function [12[ 
computed using periodic boundary conditions (k = 0). 

There are three different independent parameters (t 
sets the unit of energy) that we examine when analyz- 
ing the spectral densitiy: the interaction U, the disor- 
der A and the doping p. We start by considering the 
dependence on A for the half-filled case and constant 
[7 = 8. On increasing A, for A < U, the Mott gap 
present at A = is gradually filled in, and the evolu- 
tion of the spectral density follows the expected trend 
(see left panel of Fig. Q}. For A > U one might expect 
the residual dip to disappear completely, at least at large 
enough disorder, but this is not what is found numeri- 
cally. The right panel of Fig. Q] shows the behavior for 
A = 10, 12, 14, 16: Increasing the disorder above A = U 
leaves a residual pseudogap independent of the disorder 
strength. There is a sharpness in the behavior of A(u>) at 
small cu which suggests that the anomaly could be non- 
analytic at T = 0. However, finite size rounding prevent 
a precise characterization of this feature. 

The appearance of a pseudogap is remarkably different 
from what is found using dynamical mean field theory 
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FIG. 3: Scaling of the pseudogap as a function of interaction 
and hopping. The key observation is that the pseudogap is 
unchanged when U increases from U = 8 to U — 16 (left 
panel), and hence the antiferromagnetic exchange J is halved, 
but has a width which is proportional to t (right panel). The 
shifting in the third data set is done so as to make the two 
minima, for U = 8 and U = 16, coincide. In the last data set 
the energy axis was rescaled by a factor of two. 



(DMFT) in a similar range of parameters [13j for the same 
Hamiltonian. DMFT predicts, at least in the Hubbard- 
I approximation [l3T| . that A(lo) evolves smoothly as the 
energy crosses the chemical potential, with no pseudogap. 
This indicates that the anomaly is likely to be determined 
by non-local, short-ranged correlations and connected to 
reduced dimensionality. 

Let us now move to the incommensurate filling p = 0.6, 
where there is no Mott gap. A(lo) is plotted in Fig. [21 A 
pseudogap is now evident at all disorder strengths. For 
A > U it behaves as for the p = 1 case and saturates 
for strong enough disorder. In the A < U regime the 
pseudogap gets deeper with increasing disorder but its 
width remains largely unchanged. In contrast, the depth 
of the pseudogap which develops out of the Mott phase 
at p = 1 decreases as A increases (l7j. Such a difference is 
not surprising since, at p = 1, there is a cross-over from 
a U > A regime, dominated by the Mott-gap scale to 
the A > U regime characterized by the pseudogap scale. 
On the other hand, away from half filling, the pseudogap 
scale is always the leading instability and continues to 
persist in the highly disordered regime. It remains un- 
clear why the depth of the pseudogap saturates instead of 
moving monotonically to as A increases. Certainly the 
behavior of the anomaly away from half-filling makes evi- 
dent that the naive picture in which the effect of disorder 
is only to smear the spectral density is inappropriate. 

We have carefully studied the doping dependence of 
the pseudogap, especially since in the clean model it is 
present only for small doping. In contrast, for A = 10 



and (7 = 8 the spectral density shows a nearly universal 
pseudogap: A(uj) for p = 0.6 — 1.0 coincide over the entire 
pseudogap region. This feature persists in the A < U 
regime as soon as the system is sufficiently away from 
the Mott-gap dominated regime. Typically, at p = 0.9, 
a distinct dip at uo — is seen at all disorder strengths 
considered in this work. At densities p < 0.6 and U < 5 
the sign problem in the DQMC method is mild and does 
not prohibit obtaining accurate spectral functions [l8l ] . 
The agreement between the two numerical techniques at 
p = 0.6 and A = 8 (see right panel of Fig(2J) is excellent 
and gives a firmer basis to our speculations on the irrel- 
evance of the cluster size that was exactly diagonalized. 
DQMC results from p = 0.6 to p = 0.4 also confirms the 
stability of the pseudogap under doping. That strong 
disorder can stabilize the pseudogap was also observed 
experimentally Q in a conductivity study of YBCO sam- 
ples where Cu was substituted with Zn. The T — p phase 
diagram of this heavily disordered system shows a flat 
pseudogap crossover line in agreement with our finding 
of a doping-independent pseudogap energy scale. 

That variation in A or p both leave the anomaly un- 
changed can be understood as follows. Consider, for in- 
stance, the half-filled case with t = and A > U in which 
a finite fraction of the sites remains empty. When t ^ 0, 
electrons lower their energy by dclocalizing on neighbor- 
ing empty sites, regardless of the overall density, as the 
physics is primarily local. The effect of doping does not 
alter this situation but merely shifts the "action" on the 
new set of sites lying close to p, . Although the details of 
the effect of t are ultimately responsible for the forma- 
tion of the pseudogap, one can see that if the pseudogap 
forms and is stable upon variation in disorder then it fol- 
lows that the same anomaly forms and is stable upon 
variation of doping. 

The distinctive feature of the strongly disordered 
Anderson-Hubbard Hamiltonian is that particle localiza- 
tion occurs independently of the doping level. In the 
non-disordered case, localization is induced by correla- 
tion and gradually disappears as the density moves away 
from half-filling. Despite this different doping depen- 
dence, there are similarities with the ordered case that 
point to a common localization-related origin of the pseu- 
dogap. We found that the pseudogap scale is independent 
of U for very large U with A(ui = p) > at any finite 
doping. To make this point more quantitative we show 
the U = 8 and U = 16 pscudogaps in the left panel of 
Fig. [3] at the common disorder strength A = 8. One can 
see that, apart from a deepening of A(lo) with increasing 
U , the pseudogap scale is essentially unchanged. Inde- 
pendence of the pseudogap from U and A in the strongly 
disordered and strongly interacting regime leaves t as the 
only energy scale relevant to the pseudogap phenomenon. 
Indeed, halving t (right panel of Fig[3|) induces an almost 
linear reduction on the energy scale characterizing the 
anomaly. This suggests a kinetic mechanism for pscudo- 
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FIG. 4: The momentum resolved spectral function A^{uj) 
exhibits a pseudogap as the Fermi surface is crossed. Here 
U = 8, A = 8, and p = 0.8. 



gap formation and it certainly rules out super-exchange 
since the magnetic energy scale J tends to in this limit. 
Spin, however, has to play a crucial role since, obviously, 
spinlcss fcrmions interacting through a local U would 
show a smooth A{oj) with no suppression at the Fermi 
surface. 

An interesting feature of the clean Hubbard model and 
of experiments on cuprates is the momentum dependence 
of the pseudogap. In the hole doped case for example, 
sharp excitations survive as one crosses the Fermi energy 
along the nodal direction whereas a pseudogap develops 
along the anti-nodal one. The behavior of the averaged 
Ak (w) in the strongly disordered regime (A = 8) is given 
in Fig. [4j As k cuts through the Fermi surface either 
along the anti-nodal ((0,0) to (it, 0)) or nodal ((0,0) to 
(-7T, 7r)) lines, a depression in Ak(uj) is seen as fcp is tra- 
versed. This is therefore at odds with the results on the 
clean model and experimental data. A qualitative expla- 
nation for such a difference could reside in the different 
nature of electron localization in the two cases: while lo- 
calization is certainly isotropic in the strongly disordered 
scenario this is not necessarily so in clean materials where 
it is induced by dynamical inhomogencitics. 

In conclusion, we have shown the formation of a robust 
pseudogap in systems with strong repulsive local inter- 
action and strong disorder. This parameter regime lies 
outside the range of applicability of the perturbative AA 
calculation. It also lies outside the Coulomb gap scenario 
since the potential is local, and there is no Coulomb gap 
for the on-site Hubbard interaction when the itinerancy 
of the electrons is switched off. 

Although our exact Lanczos analysis is carried out 
only on clusters small enough to be exactly diagonalized, 
the phenomenon is also present on much larger clusters 
treated with the exact DQMC method. The pseudogap 
energy scale is set by t, a result shared with other recent 



numerical studies of non-disordered, strongly interacting 
systems 0, In particular, it persists even when J — > 0, 
and so does not appear to be linked to antiferromagnetic 
fluctuations. Finally, it is suggestive that recent experi- 
ments [f| observed a doping independent anomalous be- 
havior in the conductivity of highly disordered cuprates 
consistent with what is reported here. 

In the uniform Hubbard model, and the cuprate mate- 
rials which it may describe, spatial inhomogeneities arise 
spontaneously, without any explicit symmetry-breaking 
in the Hamiltonian itself 15j, Il6|. Likewise, the pseudo- 
gap is a feature of the model [JUL 0] and the materials. 
Both phenomena disappear with doping. Disorder how- 
ever induces i)the stabilization of the pseudogap over a 
much larger range of densities; ii) its independence on 
the momentum. These facts, together with the persis- 
tence of the pseudogap at large values of U in both dis- 
ordered and uniform models, are compatible with a sce- 
nario where the pseudogap in the strong-coupling regime 
is intrinsically linked to electron localization and driven 
by a kinetic mechanism. Such a picture, which invokes 
spatial inhomogencitics, thus connects the appearance 
of the pseudogap in the two seemingly rather different 
contexts of metal-insulator transitions driven by the in- 
terplay of randomness and correlation, and the cuprate 
superconductors. 
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